# -*- coding: utf-8 -*-
from pylab import *
import scipy.optimize as sp
from pitubalib import *

rc('font',family='serif')
rc('font',serif='Times, Times New Roman, Bitstream Vera Serif')
rc('text', usetex=False)
rc('font', size=10)
rc('xtick', labelsize=12)
rc('ytick', labelsize=12)
rc('legend', fontsize=10)
rc('figure',figsize=[5,4])


#---------------------------------------------------------------------------
# PVT FROM SERGI RESERVOIR SAMPLE
#---------------------------------------------------------------------------

## PVT 1
fluido1=classPVT()
fluido1.pe=array([200.,170.,140.,110.,80.,60.,48.,35.,25.,15.,1.033])
fluido1.Boe=array([1.165,1.149,1.133,1.116,1.100,1.089,1.083,1.076,1.069,1.062,1.0000])
fluido1.Rse=array([66.9,57.8,48.6,39.4,30.2,24.1,20.4,16.4,13.1,9.4,0.0000])
fluido1.Bge=array([0.0047,0.055,0.0067,0.0085,0.0117,0.0156,0.0195,0.0310,0.0435,0.0720,1.0000])
fluido1.Bwe=1.0*ones(len(fluido1.pe))
fluido1.viscoe=array([5.81000,7.23000,8.05000,9.05000,10.23000,11.58000,12.17000,12.80000,14.12000,14.83000,15.57000,16.33000,17.12000,17.96000,18])
fluido1.viscwe=1.0*ones(len(fluido1.pe))
fluido1.co=93e-6
fluido1.Pb=48.
fluido1.T=273.+57.

## PVT 2
fluido2=classPVT()
fluido2.pe=array([168.1700,140.1500,112.1200,84.0900,56.0600,50.4500,44.8500,39.2400,33.6300,28.0300,24.200,16.8200,11.2100,5.6100,1.0330])
fluido2.Boe=array([1.1846,1.1616,1.1385,1.1155,1.0925,1.0874,1.0828,1.0786,1.0746,1.0705,1.0660,1.0608,1.0543,1.0457,1.0000])
fluido2.Rse=array([57.3700,48.0300,38.6867,29.3433,20.0000,18.1300,16.2300,14.3500,12.4700,10.5900,8.7200,6.8400,4.9200,2.9100,0.0000])
fluido2.Bge=array([0.0015,0.0022,0.0041,0.0076,0.0193,0.0216,0.0245,0.0283,0.0333,0.0404,0.0509,0.0686,0.1039,0.2097,1.0000])
fluido2.Bwe=1.0*ones(len(fluido2.pe))
fluido2.viscoe=array([5.81000,7.23000,8.05000,9.05000,10.23000,11.58000,12.17000,12.80000,14.12000,14.83000,15.57000,16.33000,17.12000,17.96000,18])
fluido2.viscwe=1.0*ones(len(fluido2.pe))
fluido2.co=93e-6
fluido2.Pb=48.0
fluido2.T=273.+57.

fluido=fluido1


#---------------------------------------------------------------------------
# COMPUTING EQUATION COEFFICIENTS
# Np*(Bo2+(Rp-Rso2)*Bg2)+Wp*Bw2=N*(Eo+m*Eg+Efw)+(Wi+We)*Bw2+Gi*Bg2
#---------------------------------------------------------------------------
Eq=classEBM()
prod=genfromtxt('prod_SG1234.txt',delimiter='\t',skiprows=1)
Eq.Np=prod[:,0]*1e3
Eq.Gp=0.1*prod[:,1]*1e3
Eq.Gi=0.*prod[:,1]
Eq.Wp=prod[:,2]*1e3
Eq.Wi=prod[:,3]*1e3
Eq.dt=prod[:,4]

#
#Eq.Np[-1]=Eq.Np[-2]+70.*365*5
#Eq.Gp[-1]=Eq.Gp[-2]+70.*365*5*10.
#Eq.Wp[-1]=Eq.Wp[-2]+224.*365*5
#Eq.Wi[-1]=Eq.Wi[-2]+484.*365*5*1.5
#

Eq.We=zeros(len(prod[:,0]))
Eq.N=5.23e6

SG=classRes()
# DEFINING DATUM=-790M (AVERAGE RESERVOIR DEPTH)
# Pasta do poco CER-4: pe=91.3 kgf/cm2 @ MD=940m, COTA=-835.4
SG.p0= 100. #91.3+0.08*(-835.4+800.)
SG.Pb0=fluido.Pb

SG.phi=0.2
SG.Sw0=0.25
SG.Swi=0.25
SG.Soi=0.75
SG.Vp0=Eq.N*fluido.Bo(SG.p0,SG.Pb0)/(1.-SG.Sw0-SG.Sg0)


# Petrophysics
SG.m=0.
SG.cr=200e-6 #3.5e-6 * 14.7/1.033#5e-5
petroPhys=classPetrophys()
Swe=array([0.13,0.23,0.34,0.44,0.54,0.65,0.75])
petroPhys.Swe=Swe

# Aquifer
aquifero=classAquifero()
aquifero.WW=10*Eq.N # 400e6*0.27 # 35e6
aquifero.p0=SG.p0
aquifero.cr=SG.cr
aquifero.cw=1e-5
aquifero.pvt=fluido
aquifero.k=0.35
aquifero.L=1000.
aquifero.model='FetPcte'

# TESTES DE PRESSAO

# DATUM=-850m
SG.teste.p=loadtxt('testes_pressao.dat',skiprows=1)[:,1]
SG.teste.Np=loadtxt('testes_pressao.dat',skiprows=1)[:,0]*1e6
SG.teste.Gp=interp(SG.teste.Np,Eq.Np,Eq.Gp)
SG.teste.Wp=interp(SG.teste.Np,Eq.Np,Eq.Wp)
SG.teste.Wi=interp(SG.teste.Np,Eq.Np,Eq.Wi)

#
PAq=loadtxt('testes_pressao_aquifero.dat',skiprows=1)

Eq.res=SG
Eq.aquifero=aquifero
Eq.pvt=fluido
Eq.petroPhys=petroPhys

#-----------------------------------------------------------------------
# RESULTADOS DO HISTORICO
#-----------------------------------------------------------------------
# COMPUTE GAS PRODUCTION FROM EXP PRESSURE
#Gp=loadtxt('Gp.dat')
Gp=Eq.calcGp()
Eq.Gp=interp(Eq.Np,Gp[:,0],Gp[:,1])
savetxt('Gp.dat',Gp)

# COMPUTE MATERIAL BALANCE
out=classHist()
out=Eq.runHist(1.,SG.p0)

fid=open('out_prod.txt','w')
for i in range(len(prod)):
	fid.write('\n'+str(Eq.Np[i])+'\t'+str(Eq.Gp[i])+'\t'+str(Eq.Wp[i])+'\t'+str(Eq.Wi[i])+'\t'+str(out.We[i]))
fid.close()

xy=loadtxt('Np_press.dat')

figure()
plot(Eq.Np*1e-6,out.p,'k-')
plot(Eq.Np*1e-6,out.Pb,'k--')
plot(Eq.Np*1e-6,out.Pa,'k:')
plot(PAq[:,0]*1e-6,PAq[:,1],'ko')
plot(xy[:,0],xy[:,1],'ks')
#errorbar(SG.teste.Np*1e-6,SG.teste.p,fmt='ro',yerr=60./14.23)
title('Balan\c{c}o de Materiais / SG1-4')
legend(('$P_{res}\ Calc$','$P_{sat}\ Calc$','$P_{aq}\ Calc$','$P_{aq}\ Exp$', '$P_{res}\ Exp$'),'center right').draw_frame(False)
xlabel('Np ($10^6\ std\ m^3$)')
ylabel('$P_{res} (kgf/cm2)$')
ylim([0,100])
grid(1)
subplots_adjust(bottom=0.14, left=0.14)
savefig('matbal_SG1234.png',dpi=600)

figure()
plot(Eq.Np*1e-6,gradient(Eq.Gp)/gradient(Eq.Np),'rs-')
plot(Eq.Np*1e-6,zeros(len(Eq.Np)),'k:')
title('Balan\c{c}o de Materiais / SG1-4')
xlabel('Np (1e6 std m3)')
ylabel('GOR')
grid(1)

figure()
plot(Eq.Np*1e-6,Eq.Gp*1e-6,'rs-')
title('Balan\c{c}o de Materiais / SG1-4')
xlabel('Np (1e6 std m3)')
ylabel('Gp (1e6 std m3)')
grid(1)

figure()
plot(Eq.Np*1e-6,array(Eq.Wi)*1e-6,'cx-')
plot(Eq.Np*1e-6,array(out.We)*1e-6,'bx-')
plot(Eq.Np*1e-6,(-array(Eq.Wp)+array(out.We)+array(Eq.Wi))*1e-6,'rx-')
#title('Balan\c{c}o de Materiais / SG1-4.')
legend(('Wi','We','Wi+We-Wp'),'upper left')
xlabel('Np (1e6 std m3)')
ylabel('Water Volume (1e6 std m3)')
grid(1)
#ylim([0,100])

figure()
plot(Eq.Np*1e-6,out.Sw,'ks-')
plot(Eq.Np*1e-6,out.So,'kx--')
plot(Eq.Np*1e-6,out.Sg,'ko-')
plot(Eq.Np*1e-6,array(out.Sw)+array(out.So)+array(out.Sg),'-k')
xlabel('Np ($10^6\ std\ m^3$)')
legend(('Sw','So','Sg','Sw+So+Sg'),'upper right').draw_frame(False)
ylim([0,1])
grid(1)
#title('Balan\c{c} de Materiais / SG1-4')
ylabel('Saturacao de Fluidos')
savefig('matbal_saturacao.png',dpi=600)

figure()
#plot(Eq.Wp/(Eq.Np+Eq.Wp),gradient(Eq.Wp)/(gradient(Eq.Wp)+gradient(Eq.Np)),'-sk')
plot(1.-array(out.So),out.fw,'-sk')
ylabel('Fractional Flow')
xlabel('Water Saturation (Sw)')
ylim([0,1])
#legend(('Sw','So','Sg','Sw+So+Sg'),'upper right')
grid(1)
title('Fractional Flow')
#-----------------------------------------------------------------------
print 'OGIP [10^6 std m3] = ', Eq.N*fluido.Rs(SG.p0,SG.Pb0)/1e6
print 'VGIP [10^6 std m3] = ',SG.Vp(out.p[-1])*out.Sg[-1]/fluido.Bg(out.p[-1])/1e6 +(Eq.N-Eq.Np[-1])*fluido.Rs(out.Pb[-1],out.p[-1])/1e6
print 'Rs(Pb) = ', fluido1.Rs(out.Pb[-1],out.Pb[-1])
print 'Rs(P) = ', fluido1.Rs(out.p[-1],out.Pb[-1])
print 'Gas Vol RC [10^6 m3] = ',SG.Vp0*out.Sg[-1]/1e6
print 'Gas Vol ST [10^6 m3] = ',SG.Vp0*out.Sg[-1]/1e6/fluido.Bg(out.p[-1])
print 'Porous Volume [10^6 m3] = ',SG.Vp(out.p[-1])/1e6,SG.Vp0/1e6
print 'Bg = ',fluido.Bg(out.p[-1])
print 'Sg = ',out.Sg[-1]
show()
